FFT86 Fast Fourier Transform Library ==================================== Manual version 1.1A for FFT86 version 1.1 ----------------------------------------- Copyright © 1994 Vectorplex Systems. All rights reserved. FFT86 is a trademark of Vectorplex Systems. Table of Contents ----------------- Description Registering as an FFT86 User List of Files in Package Library Files FFT86E.LIB FFT867.LIB Assembler Source Code Files Function Descriptions Naming Conventions for Functions List of Functions C Header File FFT86.H FFT Error Return Code In-Place Processing of FFT Data Arrays FFT Data Scaling Difference between FFT86 Version 1.1 and 1.00 Complex-to Complex FFTs Complex Array Format cfft2f cfft2d cfft4f cfft4d Spectral Frequencies Programming Example Real-to-Complex FFTs Real Array Format RFFT Complex Format rfft2f rfft2d rfft4f rfft4d Spectral Frequencies Programming Example Window Weighting Functions Format of Array from Window Weighting Functions Arguments for Window Functions wbartltf wbartltd whammf whammd whannf whannd wkaiserf wkaiserd wparzenf wparzend Programming Example Programming Tips Get a Math Coprocessor Use Radix-2 FFTs Use Real-to-Complex FFTs Call Window Functions Only Once Float versus Double Maximum Array Length Using Maximum Length Arrays Using 8088 Class Processors Compatibility with VECTOR87 Technical Details FFT Execution Times External References Segment Address Arithmetic Modification of Code Segment Contents Windows 3.x and OS/2 Compatibility Protected Mode Operation Where is the 80386 Library? License Warranty Description ----------- FFT86 is an object-code library of 18 functions to perform the Fast Fourier Transform (FFT) and related operations. The library functions may be called from programs compiled with Microsoft C/C++. Borland C++ and other compilers that use the Microsoft C call protocol can also be used. All functions are supplied in versions that are compatible with C "float" and "double" data types. ALL functions are written in assembler language to increase execution speed and accuracy, and to reduce object code size. A complex-to-complex Cooley-Tukey FFT is provided, as well as a real-to-complex FFT. Forward and inverse transforms can be done. There are radix-2 and radix-4 versions of the FFT's. Also there is a set of functions for "weighting", or "window tapering" of time-series data before an FFT, including Hann, Hamming, Kaiser, Parzen, and Bartlett. The main processor must be an Intel 80286, 80386, or 80486, or a compatible device. If you need to use an 8088 class processor, refer to the "Programming Tips" section. An Intel compatible math coprocessor, such as the 80287, 80387, or 80487SX is optional, but highly recommended, since it will increase the speed of the FFT by a factor of 10 or more. Note that the Intel 80486DX and Pentium chips have a built-in math coprocessor. A real-to-complex FFT with an array of 1024 floats as input takes about 27 milliseconds on a 25 MHz 80486DX computer. The FFT86 library is "shareware", which means you may evaluate it for 30 days, however a license fee of $100 must be paid after this period. There are no royalty fees on applications (see the "License" section and the file "order.txt" for details). The package complies with guidelines of the ASP (Association of Shareware Professionals). For additional information, please telephone (604) 858-0832 between 9 AM and 8 PM Pacific time, send e-mail to CompuServe 70322,3532 or write to: Barrie Walker Vectorplex Systems 7572 Sapphire Drive Sardis, B.C. V2R 3A7 Canada Registering as an FFT86 User ---------------------------- We hope you find that FFT86 is useful in your work. You are allowed to evaluate the software for 30 days. Use of the software after the evaluation period requires that you register. If you register and pay the $100 license fee, you will receive: - a printed manual - a diskette containing the latest software version - full assembler source code - C source code for radix-2 FFT's such as rfft2f.c - C source code for the data windowing functions such as whannf.c - technical assistance if required Similar commercial software sells for up to $250 without source code. Although this is a small package, it took several months of work to get it fully optimized, debugged, and documented. License fees help pay for the original programming as well as updates, technical support, and development of related new products. List of Files in Package ------------------------ manual.txt this instruction manual (hardcopy is available) order.txt registration and order form readme.txt contains information not in manual fft86.h C header file for library functions fft86e.lib object library, math coprocessor optional fft867.lib object library, math coprocessor required fft86t3.c sample C program to call cfft2d fft86t4.c sample C program to call rfft2f fft86t5.c sample C program to call wkaiserf Registered users receive the following additional files on diskette: *.asm assembler source files compatible with Microsoft MASM 5.10a rfft2f.c C version of cfft2.asm, real-to-complex FFT vdsl03f.c C version of vdsl03.asm, real transform pass for FFT cfft2f.c C version of cfft2.asm, complex-to-complex FFT rfft2d.c C version of cfft2.asm, real-to-complex FFT vdsl03d.c C version of vdsl03.asm, real transform pass for FFT cfft2d.c C version of cfft2.asm, complex-to-complex FFT w*.c C versions of windowing functions fft86t2.c C program to test all functions time1.obj timer for fft86t2.c masmrun.bat assembles source code, makes .LIB files manual.doc manual in Microsoft Word for Windows format. megaflop.exe DOS program to measure floating point speed megaflop.txt manual for megaflop.exe program Library Files ============= Two versions of the FFT86 object code library are included. FFT86E.LIB is a coprocessor-optional library and FFT867.LIB is a coprocessor-required library. Select the library depending on the availability of a math coprocessor and the type of main C library that you plan to use. The library file is read by the LINK program that is part of your C compiler system. FFT86E.LIB ---------- FFT86E.LIB is the "emulator/coprocessor" version of the library. It will use a coprocessor if available; otherwise floating point instructions will be emulated by software. Use FFT86E if you want to run your program on a computer which may or may not have a coprocessor. Object code in this library was assembled using the Microsoft Macro Assembler "/E" option. FFT86E.LIB should be used in conjunction with a compatible emulator version of the main C library, which usually has a name ending with "E", like SLIBCE, CLIBCE, MLIBCE, or LLIBCE, depending on how you installed your C compiler. FFT867.LIB ---------- FFT867.LIB is the "coprocessor required" version of the library. It will work only on systems which have a math coprocessor. The code in FFT867 is smaller and floating point instructions run about 6% faster than the code in FFT86E.LIB. FFT867.LIB should be used along with a compatible coprocessor required version of the main C library, which usually has a name ending with "7", like SLIBC7, CLIBC7, MLIBC7, or LLIBC7. Assembler Source Code Files --------------------------- Note that a single assembler source code file is used to generate both the float and double version of each function. For example, rfft2.asm is used to generate rfft2f.obj and rfft2d.obj. This is done through the use of a macro and the DOS command line used to invoke the assembler. See the files masmrun.bat and macrobw2.asm for details. Function Descriptions ===================== Naming Conventions for Functions -------------------------------- The last character in all function names is either "f" or "d", to identify the float or double version. The window functions, used for Hann, Kaiser, and similar data tapering operations, have names starting with the character "w". For the FFT functions, the first 4 characters in the function name are "rfft" or "cfft", to identify the real-to-complex or complex- to-complex transform. The second last character in the FFT function name is "2" or "4", to indicate a radix-2 or radix-4 transform. List of Functions ----------------- cfft2f - complex-to-complex radix-2 FFT, float array cfft2d - complex-to-complex radix-2 FFT, double array rfft2f - real-to-complex radix-2 FFT, float array rfft2d - real-to-complex radix-2 FFT, double array cfft4f - complex-to-complex radix-4 FFT, float array cfft4d - complex-to-complex radix-4 FFT, double array rfft4f - real-to-complex radix-4 FFT, float array rfft4d - real-to-complex radix-4 FFT, double array whammf - make Hamming window, float array whammd - make Hamming window, double array whannf - make Hann window, float array whannd - make Hann window, double array wbartltf - make Bartlett window, float array wbartltd - make Bartlett window, double array wkaiserf - make Kaiser window, float array wkaiserd - make Kaiser window, double array wparzenf - make Parzen window, float array wparzend - make Parzen window, double array C Header File FFT86.H --------------------- FFT86.H is the C header file containing function prototypes for the library, and it should be "#included" in your C code. The library is compatible with the large Microsoft C memory model, however you can call the functions using any other memory model if you include FFT86.H, since the C compiler will insure that appropriate calls are done. It is suggested that you put a copy of FFT86.H in the directory where your standard C compiler library header files are located; the directory name could be C:\C600\INCLUDE depending on your compiler setup. Here is a listing of FFT86.H: // fft86.h - C function declaration header file for FFT86 library int far cfft2f (float far *pfarray, int n, int idir); int far cfft2d (double far *pdarray, int n, int idir); int far rfft2f (float far *pfarray, int n, int idir); int far rfft2d (double far *pdarray, int n, int idir); int far cfft4f (float far *pfarray, int n, int idir); int far cfft4d (double far *pdarray, int n, int idir); int far rfft4f (float far *pfarray, int n, int idir); int far rfft4d (double far *pdarray, int n, int idir); void far whammf (float far *pfarray, int n); void far whammd (double far *pdarray, int n); void far whannf (float far *pfarray, int n); void far whannd (double far *pdarray, int n); void far wbartltf (float far *pfarray, int n); void far wbartltd (double far *pdarray, int n); void far wparzenf (float far *pfarray, int n); void far wparzend (double far *pdarray, int n); void far wkaiserf (float far *pfarray, int n, float falpha); void far wkaiserd (double far *pdarray, int n, float falpha); FFT Error Return Code --------------------- All the FFT functions return an int error code, as follows: 0 - normal operation 2 - illegal FFT length 3 - array address offset nonzero for maximum length FFT In-Place Processing of FFT Data Arrays -------------------------------------- All the FFT functions do "in-place" processing. In other words, the output data array will be written over the input data array. The input data is destroyed. The output array is always exactly the same size as the input array. FFT Data Scaling ---------------- The FFT functions apply a scale factor so that data is scaled as indicated in standard formulas for the forward and inverse Discrete Fourier Transforms (DFT). The scale factor is 1/n for inverse complex-to-complex transforms and 1/(4*n) for inverse real-to- complex transforms. If you do a forward transform on a data array, followed by an inverse transform on the output from the forward transform, then the original data will be restored. Difference between FFT86 version 1.1 and 1.00 ---------------------------------------------- Version 1.00 required a password argument on the calls to all of the FFT functions. This argument is no longer needed and should not be used with version 1.1. Complex-to Complex FFTs ======================= The complex-to-complex FFT functions are cfft2f, cfft2d, cfft4f, and cfft4d. They transform n complex numbers in the input array to n complex numbers in the output array. Complex Array Format -------------------- The complex-to-complex FFTs use the same array of complex numbers for input and output. A complex array is an array of floats or doubles, as required by the function. Pairs of adjacent floats or doubles represent complex numbers, with the first float or double containing the real part, and the second float or double containing the imaginary part. For a complex-to-complex transform of length n, there must be n complex numbers in the array, represented by 2*n floats or 2*n doubles as required by the function. cfft2f ------ cfft2f is a complex-to-complex radix-2 FFT that works with an array of float data. C prototype: int far cfft2f (float far *pfarray, int n, int idir); pfarray is a far pointer to the array of float data which is to be transformed. The output data will overlay the input data. n is the transform length, which is the number of complex numbers in the input and output arrays. n must be a power of 2 between 4 and 8192, e.g. 4,8,16,...2048,4096,8192. If n is 8192 the address offset of pfarray must be zero (see "Programming Tips" section). There must be 2*n floats in the input array, since each complex number is represented by 2 floats. idir is the transform direction, which is 1 for a forward transform, and -1 for an inverse transform. A forward transform is used to transform data sampled in time to data sampled in frequency. An inverse transform is used to transform data sampled in frequency to data sampled in time. cfft2d ------ cfft2d is a complex-to-complex radix-2 FFT that works on an array of double data. cfft2d is the same as cfft2f, except an array of doubles is transformed, and the maximum transform length is reduced to 4096. C prototype: int far cfft2d (double far *pdarray, int n, int idir); pdarray is a far pointer to the array of double data which is to be transformed. The output data will overlay the input data. n is the transform length, which is the number of complex numbers in the input and output array. n must be a power of 2 between 4 and 4096, e.g. 4,8,16,...1024,2048,4096. If n=4096 the address offset of pfarray must be zero (see "Programming Tips" section). There must be 2*n doubles in the input array, since each complex number is represented by 2 doubles. idir is the transform direction, which is 1 for a forward transform, and -1 for an inverse transform. A forward transform is used to transform data sampled in time to data sampled in frequency. An inverse transform is used to transform data sampled in frequency to data sampled in time. cfft4f ------ cfft4f is the same as cfft2f, except a radix-4 algorithm is used. See the description of cfft2f for more information. cfft4d ------ cfft4d is the same as cfft2d, except a radix-4 algorithm is used. See the description of cfft2d for more information. Spectral Frequencies -------------------- For a complex-to-complex FFT, the frequencies of the spectral data may be calculated by C code similar to the fragment below: int k, n=8; /* set FFT length n */ float dt=2.0; /* set time sample interval dt=2 seconds*/ float *pfreq; /* pointer to array of frequencies in Hertz */ pfreq = (float *) calloc (n, sizeof(float)); /*get frequency array*/ for (k=0; k #include #include double darray[8] = {1.,0., 0.,0., -1.,0., 0.,0.}; /*lowest sine wave*/ int n = 4, iforward = 1, iinverse = -1, iretcode; void printcmplx (char * ctitle, int n, double far * pdarray, int icode) { int i; printf ("n=%d, icode=%d. %s \n", n, icode, ctitle); for (i=0; i<(n*2); ++i) {printf ("%5.2f ", *(pdarray+i) );} printf ("\n"); return; } void main () { printf("Test starting\n"); printcmplx ("Input to cfft2d forward:",n,darray,0); iretcode = cfft2d (darray, n, iforward); printcmplx ("Output from cfft2d forward:",n,darray,iretcode); iretcode = cfft2d (darray, n, iinverse); printcmplx ("Output from cfft2d inverse:",n,darray,iretcode); printf ("Test ending\n"); return; } The following output was produced by the C program: Test starting icode=0. Input to cfft2d forward: 1.00 0.00 0.00 0.00 -1.00 0.00 0.00 0.00 FFT86 Fast Fourier Transform Library. Version 1.00 Copyright (c) 1991 Vectorplex Data Systems Ltd. All rights reserved This unlicensed product may be evaluated for 30 days after first useage. icode=0. Output from cfft2d forward: 0.00 0.00 2.00 0.00 0.00 0.00 2.00 0.00 icode=0. Output from cfft2d inverse: 1.00 0.00 0.00 0.00 -1.00 0.00 0.00 0.00 Test ending Real-to-Complex FFTs ==================== The real-to-complex FFTs are rfft2f, rfft2d, rfft4f, and rfft4d. These functions are almost twice as fast as the corresponding complex-to-complex transforms, for the same transform length. In the forward transform direction, the real-to-complex functions transform n real values sampled in time to n/2 + 1 complex spectral values sampled in frequency stored in RFFT Complex Format. In the inverse transform direction, the functions transform n/2+1 complex spectral values sampled in frequency in RFFT Complex Format to n real values sampled in time. For inverse transforms, the functions are actually doing complex-to-real transforms. Real Array Format ----------------- The input array for a forward real-to-complex transform of length n, or the output from an inverse transform, is an array of n real numbers representing a time series. The numbers are floats or doubles, depending on the library function. RFFT Complex Format ------------------- The output array from a forward real-to-complex transform of length n is an array of n floats or doubles containing n/2 +1 complex spectral values that are compressed into the space of n/2 complex numbers. The same array format is used as input to an inverse real-to-complex FFT. The array is compressed slightly, as described below, so the input and output arrays occupy the same amount of space. This is called RFFT Complex Format. A real valued time series has no imaginary part. RFFT Complex Format makes use of three well-known properties of the Discrete Fourier Transform (DFT) of a real valued time series. First, for a transform of length n, which uses a time series containing n values, it is only necessary to calculate n/2 + 1 spectral values, ranging in frequency from zero or "D.C." to the Nyquist frequency (the Nyquist frequency is 1/(2*dt), where dt is the time series signal sampling interval in seconds). At higher frequencies, the spectral values do not need to be explicitly calculated, since they are related to the spectral values at the lower frequencies by a complex conjugate symmetry. Second, the first spectral value is at frequency zero, or D.C., and the complex spectral value has an imaginary part equal to zero. And third, the last spectral value is at the Nyquist frequency and it also has an imaginary part equal to zero. The two imaginary numbers which are known to be zero do not need to be stored in the array. RFFT Complex Format places the real part of the Nyquist spectral value in the position of the imaginary part of the D.C. spectral value. As an example, the listing below shows the format of the output array from a real-to-complex FFT, for a transform with n=8. C array index contents of array element ------- ------------------------------------------- 0 real part of first (DC) spectral value 1 real part of fifth (Nyquist) spectral value 2 real part of second spectral value 3 imaginary part of second spectral value 4 real part of third spectral value 5 imaginary part of third spectral value 6 real part of fourth spectral value 7 imaginary part of fourth spectral value In a working program, you could use the Nyquist spectral value, or you could simplify your programming by setting it to zero. Normally the spectral magnitude at the Nyquist frequency is very small since there should be no signal energy at higher frequencies, to prevent aliasing. rfft2f ------ rfft2f is a real-to-complex radix-2 FFT that works on an array of "float" data. C prototype: int far rfft2f (float far *pfarray, int n, int idir); pfarray is a far pointer to the array of float data that is to be transformed. The output data will overlay the input data. For a forward transform, the input array should contain n samples of a real time-series, and the output array will contain n/2 + 1 complex spectral values stored in RFFT Complex Format. For an inverse transform, the input array should contain n/2 + 1 complex spectral values stored in RFFT Complex Format, and the output array will contain n samples of a real time series. n is the transform length, which is the number of floats in the input and output array. n must be a power of 2 from 8 to 16384, e.g. 8,16,32,...4096,8192,16384. If n=16384 the address offset of pfarray must be zero (see "Programming Tips" section). idir is the transform direction, which is 1 for a forward transform, and -1 for an inverse transform. A forward transform is used to transform data sampled in time to data sampled in frequency. An inverse transform is used to transform data sampled in frequency to data sampled in time. rfft2d ------ rfft2d is similar to rfft2f, except an array of double data is transformed, and the maximum transform length is 8192. C prototype: int far rfft2d (double far *pdarray, int n); pdarray is a far pointer to the array of double data which is to be transformed. The output data will overlay the input data. For a forward transform, the input array should contain n samples of a real time-series, and the output array will contain n/2 + 1 complex spectral values stored in RFFT Complex Format. For an inverse transform, the input array should contain n/2 + 1 complex spectral values stored in RFFT Complex Format, and the output array will contain n samples of a real time series. n is the transform length, which is the number of doubles in the input and output arrays. n must be a power of 2 from 8 to 8192, e.g. 8,16,32,...2048,4096,8192. If n=8192 the address offset of pdarray must be zero (see "Programming Tips" section). idir is the transform direction, which is 1 for a forward transform, and -1 for an inverse transform. A forward transform is used to transform data sampled in time to data sampled in frequency. An inverse transform is used to transform data sampled in frequency to data sampled in time. rfft4f ------ rfft4f is the same as rfft2f, except a radix-4 algorithm is used, see the description of rfft2f for more information. rfft4d ------ rfft4d is the same as rfft2d, except a radix-4 algorithm is used, see the description of rfft2d for more information. Spectral Frequencies -------------------- For a real-to-complex FFT, the frequencies of the spectral data may be calculated by C code similar to the fragment below: int k, n=8; /* set FFT length n */ int nfreq; /*number of frequencies will be n/2+1 */ float dt=2.0; /*set time sample interval dt=2 seconds*/ float *pfreq; /* pointer to array of frequencies in Hertz */ nfreq = n/2 + 1; /*no. frequencies in RFFT Complex Format*/ pfreq = (float *) calloc (nfreq, sizeof(float)); /*frequency array*/ for (k=0; k #include #include #include float * pfarray; int n = 8, iforward = 1, iinverse = -1, iretcode; void printreal (char *ctitle, int n, float far *pfarray, int icode) { int i; printf ("n=%d, icode=%d. %s \n", n, icode, ctitle); for (i=0; i #include #include #include float *pfarray, *pfwindow; int n = 8, iforward = 1, iinverse = -1, iretcode; void printreal (char * ctitle, int n, float far * pfarray, int icode) { int i; printf ("n=%d, icode=%d. %s \n", n, icode, ctitle); for (i=0; i